Make the prediction grid

Author

Max Lindmark

Published

August 18, 2023

Intro

Make an evenly spaced UTM prediction grid with all spatially varying covariates for the diet and the biomass data

pkgs <- c("tidyverse", "tidylog", "sp", "raster", "devtools", "RCurl", "sdmTMB", "terra", "ncdf4", "chron") 

if(length(setdiff(pkgs, rownames(installed.packages()))) > 0){

    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
  
}

invisible(lapply(pkgs, library, character.only = T))

# Source code for map plots
devtools::source_url("https://raw.githubusercontent.com/maxlindmark/pred-prey-overlap/main/R/functions/map-plot.R")

# Packages not on CRAN
# devtools::install_github("seananderson/ggsidekick") # not on CRAN 
library(ggsidekick)
theme_set(theme_sleek())

# Set path
home <- here::here()

Read data and depth-raster

# Read data
d <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/catch_clean.csv") |> 
  rename(X = x, Y = y)
Rows: 9376 Columns: 26
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (5): country, haul_id, ices_rect, month_year, substrate
dbl (21): year, lat, lon, quarter, month, sub_div, oxy, temp, sal, depth, x,...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed 2 variables (X, Y)

Make the grid with depth

First make a grid for the biomass data, then subset that based on the extend of the stomach data

x <- d$X
y <- d$Y
z <- chull(x, y)

coords <- cbind(x[z], y[z])

coords <- rbind(coords, coords[1, ])

plot(coords[, 1] ~ coords[, 2]) # plot data

sp_poly <- sp::SpatialPolygons(
  list(sp::Polygons(list(sp::Polygon(coords)), ID = 1))
  )

sp_poly_df <- sp::SpatialPolygonsDataFrame(sp_poly,
                                           data = data.frame(ID = 1)
                                           )
cell_width <- 3

pred_grid <- expand.grid(
  X = seq(min(d$X), max(d$X), cell_width),
  Y = seq(min(d$Y), max(d$Y), cell_width),
  year = c(1993:2021)
  )

ggplot(pred_grid |> filter(year == 2019), aes(X, Y)) +
  geom_point(size = 0.1) +
  theme_void() +
  coord_sf()
filter: removed 882,420 rows (97%), 31,515 rows remaining

sp::coordinates(pred_grid) <- c("X", "Y")

inside <- !is.na(sp::over(pred_grid, as(sp_poly_df, "SpatialPolygons")))

pred_grid <- pred_grid[inside, ]

pred_grid <- as.data.frame(pred_grid)

ggplot(data = filter(pred_grid, year == 1999), aes(X*1000, Y*1000)) + 
  geom_point(size = 0.001, alpha = 0.5) +
  NULL
filter: removed 531,776 rows (97%), 18,992 rows remaining

plot_map +
  geom_point(data = filter(pred_grid, year == 1999), aes(X*1000, Y*1000), size = 0.001, alpha = 0.5) +
  NULL
filter: removed 531,776 rows (97%), 18,992 rows remaining
Warning: Removed 704 rows containing missing values (`geom_point()`).

# Add lat and lon
# Need to go from UTM to lat long for this one...
# https://stackoverflow.com/questions/30018098/how-to-convert-utm-coordinates-to-lat-and-long-in-r
xy <- as.matrix(pred_grid |> dplyr::select(X, Y) |> mutate(X = X*1000, Y = Y*1000))
mutate: changed 550,768 values (100%) of 'X' (0 new NA)
        changed 550,768 values (100%) of 'Y' (0 new NA)
v <- vect(xy, crs="+proj=utm +zone=33 +datum=WGS84  +units=m")
y <- project(v, "+proj=longlat +datum=WGS84")
lonlat <- geom(y)[, c("x", "y")]

pred_grid$lon <- lonlat[, 1]
pred_grid$lat <- lonlat[, 2]

ggplot(filter(pred_grid, year == 1999), aes(lon, lat)) + geom_point()
filter: removed 531,776 rows (97%), 18,992 rows remaining

# Add depth now to remove islands and remaining land
# https://gis.stackexchange.com/questions/411261/read-multiple-layers-raster-from-ncdf-file-using-terra-package
# https://emodnet.ec.europa.eu/geoviewer/
dep_raster <- terra::rast(paste0(home, "/data/Mean depth natural colour (with land).nc"))
class(dep_raster)
[1] "SpatRaster"
attr(,"package")
[1] "terra"
crs(dep_raster, proj = TRUE)
[1] "+proj=longlat +datum=WGS84 +no_defs"
plot(dep_raster)

pred_grid$depth <- terra::extract(dep_raster, pred_grid |> dplyr::select(lon, lat))$elevation

ggplot(pred_grid, aes(lon, lat, color = depth*-1)) + 
  geom_point()

pred_grid$depth <- pred_grid$depth*-1

pred_grid <- pred_grid |> drop_na(depth)
drop_na: removed 76,995 rows (14%), 473,773 rows remaining
pred_grid |> 
  filter(year == 1999) |> 
  drop_na(depth) |> 
  #mutate(water = ifelse(depth < 0.00000001, "N", "Y")) |> 
  ggplot(aes(X*1000, Y*1000, fill = depth)) + 
  geom_raster() +
  NULL
filter: removed 457,436 rows (97%), 16,337 rows remaining
drop_na: no rows removed

plot_map + 
  geom_point(data = pred_grid, aes(X*1000, Y*1000), size = 0.001) + 
  geom_sf()
Warning: Removed 20329 rows containing missing values (`geom_point()`).

plot_map + 
  geom_raster(data = filter(pred_grid, year == 1999), aes(X*1000, Y*1000, fill = depth), size = 0.001) + 
  geom_sf()
filter: removed 457,436 rows (97%), 16,337 rows remaining
Warning in geom_raster(data = filter(pred_grid, year == 1999), aes(X * 1000, :
Ignoring unknown parameters: `size`
Warning: Removed 762 rows containing missing values (`geom_raster()`).

Substrate

First add lat and lon based on X and Y (utm)

substrate <- terra::rast(paste0(home, "/data/substrate-tif/BALANCE_SEABED_SEDIMENT.tif"))

newcrs <- "+proj=longlat +datum=WGS84"

substrate_longlat = terra::project(substrate, newcrs)

|---------|---------|---------|---------|
=========================================
                                          
plot(substrate_longlat)

# Now extract the values from the saduria raster to pred_grid
pred_grid$substrate <- terra::extract(substrate_longlat, pred_grid |> dplyr::select(lon, lat))$BALANCE_SEABED_SEDIMENT

unique(pred_grid$substrate)
  [1] 3.000000 2.207703 2.000000 5.000000 2.751267 3.000000 2.411566 3.000000
  [9]       NA 5.000000 3.023438 3.832703 4.693979 4.816548 2.407304 2.675289
 [17] 5.000000 2.515299 2.104980 3.309301 2.215659 4.000000 2.418371 3.475651
 [25] 3.136963 3.113001 2.414551 4.865710 2.082968 3.555176 3.868287 3.925187
 [33] 3.036916 4.963291 2.703572 4.166504 3.295166 3.913574 2.989511 2.172908
 [41] 3.014659 3.177487 4.189453 4.185059 2.021204 3.047852 2.189942 4.947421
 [49] 2.084229 2.895818 4.771301 3.115845 2.467070 4.691040 2.189209 2.446289
 [57] 3.239676 3.565430 2.035613 4.346191 4.295410 2.221002 3.051758 2.749512
 [65] 3.481570 2.115028 2.060020 4.481808 3.367188 3.831543 1.000000 4.013517
 [73] 3.182683 2.004699 4.814331 2.605957 2.217144 1.983107 3.768313 2.091238
 [81] 3.106757 2.150879 3.954897 2.124500 3.680969 3.013173 3.387896 2.451553
 [89] 2.987671 3.004685 1.230957 3.087122 4.163172 4.331752 2.847656 3.770158
 [97] 3.125428 4.137321 3.885742 3.675119 4.481557 3.618304 3.535532 4.465820
[105] 2.394043 3.985867 3.518455 2.930088 4.570766 3.035295 4.109087 3.688411
[113] 3.070784 3.621094 2.226618 3.092285 4.952732 4.725516 3.052990 3.999512
[121] 4.006014 3.395530 3.495117 3.194580 2.990574 4.053308 3.811279 4.897705
[129] 4.969238 3.747452 4.660645 4.250865 2.874023 2.834995 2.995726 2.603516
[137] 4.352661 4.529053 2.199267 3.234183 4.856023 4.949766 3.394751 3.105818
[145] 4.119141 2.052490 3.812706 2.458618 3.877474 4.508667 4.997502 4.851562
[153] 3.038879 2.274684 2.144043 2.439031 3.721595 3.773949 3.876709 4.546921
[161] 4.142808 3.771973 3.371582 4.526794 4.744454 2.290039 2.187041 3.043457
[169] 3.958679 4.084656 2.060291 4.713074 4.014312 2.630943 4.024496 4.889182
[177] 4.974456 3.012340 4.259968 4.845192 3.074554 4.409346 3.289622 3.797133
[185] 4.741160 4.990265 4.992432 2.236874 4.020061 4.972815 3.099773 3.104312
[193] 3.131378 4.014556 4.802023 3.626464 2.661095 4.346972 3.388261 2.413086
[201] 3.486816 4.047725 4.273254 3.218506 4.880537 2.026439 3.735108 3.617798
[209] 3.750732 3.903769 3.027108 4.974083 2.000000 3.153126 4.907196 4.142573
[217] 2.127777 3.366211 3.833305 3.686989 3.969047 2.451714 3.201564 3.003133
[225] 3.956664 2.793945 3.498962 3.130957 3.051025 3.170654 3.023042 2.039551
[233] 4.997326 3.213953 2.033055 2.056885 3.334961 2.716752 2.680715 3.365484
[241] 4.719290 4.165028 3.439092 2.860311 4.233764 2.220473 3.163818 4.547607
[249] 2.650146 3.881836 2.063848 2.738892 2.771484 3.844098 4.561007 2.790447
[257] 3.166412 2.996099 3.073975 4.000000 4.284601 3.897705 2.213647 2.750740
[265] 3.372153 3.138967 4.034568 2.040667 3.311523 4.179138 2.755887 2.779238
[273] 3.002139 2.831177 4.645435 4.891513 4.455050 3.114136 2.985744 4.142618
[281] 4.969971 4.160156 2.700928 2.437500 3.261470 2.720584 2.674125 4.945322
[289] 4.977206 4.981665 2.354310 3.050293 3.362305 4.034168 3.144043 3.427612
[297] 3.115974 4.706665 4.172024 4.919174 4.779785 4.999612 3.258278 4.724852
[305] 2.403809 2.568726 4.936097 4.244507 4.219971 4.019531 4.826586 3.978732
[313] 4.504028 3.343726 4.380249 2.978496 4.912631 3.741340 3.936610 2.325732
[321] 4.247200 2.162385 2.233503 2.387978 4.661255 3.199463 3.274317 2.430664
[329] 2.529657 2.882244 3.927612 4.913910 2.105009 3.039979 3.861948 4.729654
[337] 3.608406 2.118617 2.201874 3.994564 2.977833 2.230173 3.841365 2.345091
[345] 1.502316 3.117709 3.326479 4.573112 2.669378 4.844849 3.000950 2.726405
[353] 3.340504 4.578949 4.333606 4.793458 1.626249 2.945679 3.053228 3.968295
[361] 2.776001 1.726807 3.005737 3.789629 2.304810 3.019450 4.060608 2.094116
[369] 1.936951 3.144531 3.184711 2.283091 2.409424 1.728929 4.570679 4.402838
[377] 4.916186 2.350037 4.675385 4.664307 3.292003 2.157396 3.532034 4.424982
[385] 4.074768 4.896749 4.358226 2.311279 3.726698 3.015447 3.681493 4.969622
[393] 2.488953 2.328379 3.550172 4.602893 3.802612 4.695945 4.574463 4.006093
[401] 3.893770 4.331733 3.730713 4.537598 2.354356 4.375757 4.304613 4.836199
[409] 4.999268 4.829179 3.318237 4.121094 2.585449 4.685639 4.022016 4.823242
[417] 4.928812 2.574462 2.792725 4.440264 4.267936 4.459198 4.037354 2.681641
[425] 3.636063 4.099470 3.484864 4.835382 4.030147 4.965707 3.045020 4.963444
[433] 4.151367 4.356689 4.242188 4.018444 4.086391 4.748169 4.875027 2.393384
[441] 4.611573 3.831160 3.592084 2.666638 4.245511 4.447388 3.375265 4.155386
[449] 3.733628 2.018555 4.718811 4.296753 3.305542 4.881466 4.944336 4.507011
[457] 4.219727 4.868225 2.236110 4.212225 4.623746 3.489990 4.505832
factor(sort(unique(round(pred_grid$substrate))))
[1] 1 2 3 4 5
Levels: 1 2 3 4 5
pred_grid$substrate <- round(pred_grid$substrate)

pred_grid <- pred_grid %>% mutate(substrate = ifelse(substrate == 1, "bedrock", substrate),
                                  substrate = ifelse(substrate == 2, "hard-bottom complex", substrate),
                                  substrate = ifelse(substrate == 3, "sand", substrate),
                                  substrate = ifelse(substrate == 4, "hard clay", substrate),
                                  substrate = ifelse(substrate == 5, "mud", substrate))
mutate: converted 'substrate' from double to character (0 new NA)
# I. Bedrock.
# II. Hard bottom complex, includes patchy hard surfaces and coarse sand (sometimes also clay) to boulders.
# III. Sand including fine to coarse sand (with gravel exposures).
# IV. Hard clay sometimes/often/possibly exposed or covered with a thin layer of sand/gravel.
# V. Mud including gyttja-clay to gyttja-silt.

# Plot
ggplot(pred_grid, aes(X, Y, fill = substrate)) +
  geom_raster() +
  coord_sf()

Add month_year variable for matching with raster layers

pred_grid_1 <- pred_grid |> mutate(quarter = 1)
mutate: new variable 'quarter' (double) with one unique value and 0% NA
pred_grid_4 <- pred_grid |> mutate(quarter = 4)
mutate: new variable 'quarter' (double) with one unique value and 0% NA
dat <- bind_rows(pred_grid_1, pred_grid_4) |> 
  mutate(month = ifelse(quarter == 1, 2, 11), # most common months
         month_year = paste(month, year, sep = "_"))
mutate: new variable 'month' (double) with 2 unique values and 0% NA
        new variable 'month_year' (character) with 58 unique values and 0% NA

Oxygen

# Downloaded from here: https://data.marine.copernicus.eu/products
# Extract raster points: https://gisday.wordpress.com/2014/03/24/extract-raster-values-from-points-using-r/comment-page-1/
# https://rpubs.com/boyerag/297592
# https://pjbartlein.github.io/REarthSysSci/netCDF.html#get-a-variable
# Open the netCDF file
ncin <- nc_open(paste0(home, "/data/NEMO/cmems_mod_bal_bgc_my_P1M-m_1691065496556.nc"))

print(ncin)
File /Users/maxlindmark/Dropbox/Max work/R/pred-prey-overlap/data/NEMO/cmems_mod_bal_bgc_my_P1M-m_1691065496556.nc (NC_FORMAT_CLASSIC):

     1 variables (excluding dimension variables):
        float o2b[lon,lat,time]   
            standard_name: mole_concentration_of_dissolved_molecular_oxygen_in_sea_water
            long_name: Dissolved_Oxygen_Concentration
            units: mmol m-3
            _FillValue: -999
            missing_value: -999
            cell_methods: time: mean
            unit_long: millimole O2 per unit volume
            interval_write: 1 d
            interval_operation: 90 s
            online_operation: average
            _ChunkSizes: 1
             _ChunkSizes: 774
             _ChunkSizes: 763

     3 dimensions:
        time  Size:347 
            standard_name: time
            long_name: time
            bounds: time_bnds
            units: days since 1900-01-01 00:00:00
            calendar: standard
            axis: T
            _ChunkSizes: 512
            _CoordinateAxisType: Time
            valid_min: 34013
            valid_max: 44544.5
        lat  Size:498 
            standard_name: latitude
            long_name: Latitude of each location
            units: degrees_north
            axis: Y
            _CoordinateAxisType: Lat
            valid_min: 53.0082969665527
            valid_max: 61.2914962768555
        lon  Size:519 
            standard_name: longitude
            long_name: Longitude of each location
            units: degrees_east
            axis: X
            _CoordinateAxisType: Lon
            valid_min: 10.5414848327637
            valid_max: 24.9298248291016

    22 global attributes:
        CDI: Climate Data Interface version 1.9.9rc1 (https://mpimet.mpg.de/cdi)
        Conventions: CF-1.0
        history: Tue Mar 14 13:31:41 2023: cdo -z zip_1 --sortname copy /net/isilon/ifs/arch/home/iri/BALMFC/Nemo4.0/BALMFCMYP2022_univariateSST//BAL-ERGOM_BGC-DailyMeans-20211201.nc tmp_20211201.nc
        source: CMEMS BAL MFC NEMO model output converted to NetCDF
        institution: Baltic MFC, PU Danish Meteorological Institute
        comment: Data on cropped native product grid. Horizontal velocities destagged
        grid_resolution: ~1 nautical mile (1min latitude; 1min40sec longitude)
        contact: servicedesk.cmems@mercator-ocean.eu
        references: https://marine.copernicus.eu/
        file_quality_index: 1
        compression: yes
        stop_date: 2021-12-01 12:00:00
        creation_date: 2023-03-14 13:31:27
        FROM_ORIGINAL_FILE__netcdf_version_id: 4.6.1 of Oct 19 2018 20:03:53 $
        FROM_ORIGINAL_FILE__easternmost_longitude: 30.2074012756348
        FROM_ORIGINAL_FILE__northernmost_latitude: 65.8909912109375
        FROM_ORIGINAL_FILE__westernmost_longitude: 9.04154205322266
        FROM_ORIGINAL_FILE__southernmost_latitude: 53.0082969665527
        start_date: 2021-12-01 12:00:00
        CDO: Climate Data Operators version 1.9.9rc1 (https://mpimet.mpg.de/cdo)
        title: CMEMS ERGOM monthly integrated model fields
        _CoordSysBuilder: ucar.nc2.dataset.conv.CF1Convention
# Get longitude and latitude
lon <- ncvar_get(ncin,"lon")
nlon <- dim(lon)
head(lon)
[1] 10.54148 10.56926 10.59704 10.62481 10.65259 10.68037
lat <- ncvar_get(ncin,"lat")
nlat <- dim(lat)
head(lat)
[1] 53.00830 53.02496 53.04163 53.05830 53.07496 53.09163
# Get time
time <- ncvar_get(ncin,"time")
time
  [1] 34013.0 34042.5 34073.0 34103.5 34134.0 34164.5 34195.5 34226.0 34256.5
 [10] 34287.0 34317.5 34348.5 34378.0 34407.5 34438.0 34468.5 34499.0 34529.5
 [19] 34560.5 34591.0 34621.5 34652.0 34682.5 34713.5 34743.0 34772.5 34803.0
 [28] 34833.5 34864.0 34894.5 34925.5 34956.0 34986.5 35017.0 35047.5 35078.5
 [37] 35108.5 35138.5 35169.0 35199.5 35230.0 35260.5 35291.5 35322.0 35352.5
 [46] 35383.0 35413.5 35444.5 35474.0 35503.5 35534.0 35564.5 35595.0 35625.5
 [55] 35656.5 35687.0 35717.5 35748.0 35778.5 35809.5 35839.0 35868.5 35899.0
 [64] 35929.5 35960.0 35990.5 36021.5 36052.0 36082.5 36113.0 36143.5 36174.5
 [73] 36204.0 36233.5 36264.0 36294.5 36325.0 36355.5 36386.5 36417.0 36447.5
 [82] 36478.0 36508.5 36539.5 36569.5 36599.5 36630.0 36660.5 36691.0 36721.5
 [91] 36752.5 36783.0 36813.5 36844.0 36874.5 36905.5 36935.0 36964.5 36995.0
[100] 37025.5 37056.0 37086.5 37117.5 37148.0 37178.5 37209.0 37239.5 37270.5
[109] 37300.0 37329.5 37360.0 37390.5 37421.0 37451.5 37482.5 37513.0 37543.5
[118] 37574.0 37604.5 37635.5 37665.0 37694.5 37725.0 37755.5 37786.0 37816.5
[127] 37847.5 37878.0 37908.5 37939.0 37969.5 38000.5 38030.5 38060.5 38091.0
[136] 38121.5 38152.0 38182.5 38213.5 38244.0 38274.5 38305.0 38335.5 38366.5
[145] 38396.0 38425.5 38456.0 38486.5 38517.0 38547.5 38578.5 38609.0 38639.5
[154] 38670.0 38700.5 38731.5 38761.0 38790.5 38821.0 38851.5 38882.0 38912.5
[163] 38943.5 38974.0 39004.5 39035.0 39065.5 39096.5 39126.0 39155.5 39186.0
[172] 39216.5 39247.0 39277.5 39308.5 39339.0 39369.5 39400.0 39430.5 39461.5
[181] 39491.5 39521.5 39552.0 39582.5 39613.0 39643.5 39674.5 39705.0 39735.5
[190] 39766.0 39796.5 39827.5 39857.0 39886.5 39917.0 39947.5 39978.0 40008.5
[199] 40039.5 40070.0 40100.5 40131.0 40161.5 40192.5 40222.0 40251.5 40282.0
[208] 40312.5 40343.0 40373.5 40404.5 40435.0 40465.5 40496.0 40526.5 40557.5
[217] 40587.0 40616.5 40647.0 40677.5 40708.0 40738.5 40769.5 40800.0 40830.5
[226] 40861.0 40891.5 40922.5 40952.5 40982.5 41013.0 41043.5 41074.0 41104.5
[235] 41135.5 41166.0 41196.5 41227.0 41257.5 41288.5 41318.0 41347.5 41378.0
[244] 41408.5 41439.0 41469.5 41500.5 41531.0 41561.5 41592.0 41622.5 41653.5
[253] 41683.0 41712.5 41743.0 41773.5 41804.0 41834.5 41865.5 41896.0 41926.5
[262] 41957.0 41987.5 42018.5 42048.0 42077.5 42108.0 42138.5 42169.0 42199.5
[271] 42230.5 42261.0 42291.5 42322.0 42352.5 42383.5 42413.5 42443.5 42474.0
[280] 42504.5 42535.0 42565.5 42596.5 42627.0 42657.5 42688.0 42718.5 42749.5
[289] 42779.0 42808.5 42839.0 42869.5 42900.0 42930.5 42961.5 42992.0 43022.5
[298] 43053.0 43083.5 43114.5 43144.0 43173.5 43204.0 43234.5 43265.0 43295.5
[307] 43326.5 43357.0 43387.5 43418.0 43448.5 43479.5 43509.0 43538.5 43569.0
[316] 43599.5 43630.0 43660.5 43691.5 43722.0 43752.5 43783.0 43813.5 43844.5
[325] 43874.5 43904.5 43935.0 43965.5 43996.0 44026.5 44057.5 44088.0 44118.5
[334] 44149.0 44179.5 44210.5 44240.0 44269.5 44300.0 44330.5 44361.0 44391.5
[343] 44422.5 44453.0 44483.5 44514.0 44544.5
tunits <- ncatt_get(ncin,"time","units")
nt <- dim(time)
nt
[1] 347
tunits
$hasatt
[1] TRUE

$value
[1] "days since 1900-01-01 00:00:00"
# Get oxygen
dname <- "o2b"

oxy_array <- ncvar_get(ncin,dname)
dlname <- ncatt_get(ncin,dname,"long_name")
dunits <- ncatt_get(ncin,dname,"units")
fillvalue <- ncatt_get(ncin,dname,"_FillValue")
dim(oxy_array)
[1] 519 498 347
# Get global attributes
title <- ncatt_get(ncin,0,"title")
institution <- ncatt_get(ncin,0,"institution")
datasource <- ncatt_get(ncin,0,"source")
references <- ncatt_get(ncin,0,"references")
history <- ncatt_get(ncin,0,"history")
Conventions <- ncatt_get(ncin,0,"Conventions")

# Convert time: split the time units string into fields
tustr <- strsplit(tunits$value, " ")
tdstr <- strsplit(unlist(tustr)[3], "-")
tmonth <- as.integer(unlist(tdstr)[2])
tday <- as.integer(unlist(tdstr)[3])
tyear <- as.integer(unlist(tdstr)[1])

# Here I deviate from the guide a little bit. Save this info:
dates <- chron(time, origin = c(tmonth, tday, tyear))

# Crop the date variable
months <- as.numeric(substr(dates, 2, 3))
years <- as.numeric(substr(dates, 8, 9))
years <- ifelse(years > 90, 1900 + years, 2000 + years)

# Replace netCDF fill values with NA's
oxy_array[oxy_array == fillvalue$value] <- NA

dim(oxy_array)
[1] 519 498 347
str(dim(oxy_array))
 int [1:3] 519 498 347
# The third slot is the date index

# Loop through all "dates", put into a list 
dlist <- list()

for(i in 1:length(months)) {
  
  oxy_sub <- oxy_array[, , i]
    
  dlist[[i]] <- oxy_sub
  
}

# Name the list
names(dlist) <- paste(months, years, sep = "_")
str(dlist)
List of 347
 $ 2_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_1993: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_1993: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_1993: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_1994: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_1994: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_1994: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_1995: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_1995: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_1995: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_1996: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_1996: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_1996: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_1997: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_1997: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_1997: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_1998: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_1998: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_1998: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_1999: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_1999: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_1999: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_2000: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_2000: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_2000: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
  [list output truncated]
# Create data holding object
oxy_data_list <- list()

# Loop through each month_year and extract raster values for the cpue data points
for(i in unique(dat$month_year)) { # We can use q1 as looping index, doesn't matter!
  
  # Set plot limits
  ymin = 54; ymax = 58; xmin = 12; xmax = 22

  # Subset a month-year combination
  oxy_slice <- dlist[[i]]
  
  # Create raster for that year (i)
  r <- raster(t(oxy_slice), xmn = min(lon), xmx = max(lon), ymn = min(lat), ymx = max(lat),
              crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
  
  # Flip...
  r <- flip(r, direction = 'y')
  
  plot(r, main = paste(i), ylim = c(ymin, ymax), xlim = c(xmin, xmax))

  # Filter the same year (i) in the cpue data and select only coordinates
  d_slice <- dat %>% filter(month_year == i) %>% dplyr::select(lon, lat)
  
  # Make into a SpatialPoints object
  data_sp <- SpatialPoints(d_slice)
  
  # Extract raster value (oxygen)
  rasValue <- raster::extract(r, data_sp)
  
  # Now we want to plot the results of the raster extractions by plotting the cpue data points over a raster and saving it for each year.
  # Make the SpatialPoints object into a raster again (for plot)
  df <- as.data.frame(data_sp)
  
  # Add in the raster value in the df holding the coordinates for the cpue data
  d_slice$oxy <- rasValue
  
  # Add in which year
  d_slice$month_year <- i

  # Now the unit of oxygen is mmol/m3. I want it to be ml/L. The original model is in unit ml/L
  # and it's been converted by the data host. Since it was converted without accounting for
  # pressure or temperature, I can simply use the following conversion factor:
  # 1 ml/l = 103/22.391 = 44.661 μmol/l -> 1 ml/l = 0.044661 mmol/l = 44.661 mmol/m^3 -> 0.0223909 ml/l = 1mmol/m^3
  # https://ocean.ices.dk/tools/unitconversion.aspx

  d_slice$oxy <- d_slice$oxy * 0.0223909
    
  # Add each years' data in the list
  oxy_data_list[[i]] <- d_slice

}
filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

# Now create a data frame from the list of all annual values
big_dat_oxy <- dplyr::bind_rows(oxy_data_list)

Temperature

# Open the netCDF file
ncin <- nc_open(paste0(home, "/data/NEMO/cmems_mod_bal_phy_my_P1M-m_1691065935335.nc"))
                                        
print(ncin)
File /Users/maxlindmark/Dropbox/Max work/R/pred-prey-overlap/data/NEMO/cmems_mod_bal_phy_my_P1M-m_1691065935335.nc (NC_FORMAT_CLASSIC):

     1 variables (excluding dimension variables):
        float bottomT[lon,lat,time]   
            standard_name: sea_water_potential_temperature_at_sea_floor
            long_name: Sea floor potential temperature
            units: degrees_C
            _FillValue: -999
            missing_value: -999
            cell_methods: time: mean
            unit_long: degrees Celsius
            interval_write: 1 d
            interval_operation: 90 s
            online_operation: average
            _ChunkSizes: 1
             _ChunkSizes: 774
             _ChunkSizes: 763

     3 dimensions:
        time  Size:347 
            standard_name: time
            long_name: time
            bounds: time_bnds
            units: days since 1900-01-01 00:00:00
            calendar: standard
            axis: T
            _ChunkSizes: 512
            _CoordinateAxisType: Time
            valid_min: 34013
            valid_max: 44544.5
        lat  Size:498 
            standard_name: latitude
            long_name: Latitude of each location
            units: degrees_north
            axis: Y
            _CoordinateAxisType: Lat
            valid_min: 53.0082969665527
            valid_max: 61.2914962768555
        lon  Size:519 
            standard_name: longitude
            long_name: Longitude of each location
            units: degrees_east
            axis: X
            _CoordinateAxisType: Lon
            valid_min: 10.5414848327637
            valid_max: 24.9298248291016

    22 global attributes:
        CDI: Climate Data Interface version 1.9.9rc1 (https://mpimet.mpg.de/cdi)
        Conventions: CF-1.0
        source: CMEMS BAL MFC NEMO model output converted to NetCDF
        institution: Baltic MFC, PU Danish Meteorological Institute
        comment: Data on cropped native product grid. Horizontal velocities destagged
        grid_resolution: ~1 nautical mile (1min latitude; 1min40sec longitude)
        contact: servicedesk.cmems@mercator-ocean.eu
        references: https://marine.copernicus.eu/
        file_quality_index: 1
        compression: yes
        stop_date: 2021-12-01 12:00:00
        creation_date: 2023-03-14 13:30:28
        FROM_ORIGINAL_FILE__netcdf_version_id: 4.6.1 of Oct 19 2018 20:03:53 $
        FROM_ORIGINAL_FILE__easternmost_longitude: 30.2074012756348
        FROM_ORIGINAL_FILE__northernmost_latitude: 65.8909912109375
        FROM_ORIGINAL_FILE__westernmost_longitude: 9.04154205322266
        FROM_ORIGINAL_FILE__southernmost_latitude: 53.0082969665527
        start_date: 2021-12-01 12:00:00
        CDO: Climate Data Operators version 1.9.9rc1 (https://mpimet.mpg.de/cdo)
        title: CMEMS NEMO monthly integrated model fields
        _CoordSysBuilder: ucar.nc2.dataset.conv.CF1Convention
        history: Data extracted from dataset http://localhost:8080/thredds/dodsC/cmems_mod_bal_phy_my_P1M-m
# Get longitude and latitude
lon <- ncvar_get(ncin,"lon")
nlon <- dim(lon)
head(lon)
[1] 10.54148 10.56926 10.59704 10.62481 10.65259 10.68037
lat <- ncvar_get(ncin,"lat")
nlat <- dim(lat)
head(lat)
[1] 53.00830 53.02496 53.04163 53.05830 53.07496 53.09163
# Get time
time <- ncvar_get(ncin,"time")
time
  [1] 34013.0 34042.5 34073.0 34103.5 34134.0 34164.5 34195.5 34226.0 34256.5
 [10] 34287.0 34317.5 34348.5 34378.0 34407.5 34438.0 34468.5 34499.0 34529.5
 [19] 34560.5 34591.0 34621.5 34652.0 34682.5 34713.5 34743.0 34772.5 34803.0
 [28] 34833.5 34864.0 34894.5 34925.5 34956.0 34986.5 35017.0 35047.5 35078.5
 [37] 35108.5 35138.5 35169.0 35199.5 35230.0 35260.5 35291.5 35322.0 35352.5
 [46] 35383.0 35413.5 35444.5 35474.0 35503.5 35534.0 35564.5 35595.0 35625.5
 [55] 35656.5 35687.0 35717.5 35748.0 35778.5 35809.5 35839.0 35868.5 35899.0
 [64] 35929.5 35960.0 35990.5 36021.5 36052.0 36082.5 36113.0 36143.5 36174.5
 [73] 36204.0 36233.5 36264.0 36294.5 36325.0 36355.5 36386.5 36417.0 36447.5
 [82] 36478.0 36508.5 36539.5 36569.5 36599.5 36630.0 36660.5 36691.0 36721.5
 [91] 36752.5 36783.0 36813.5 36844.0 36874.5 36905.5 36935.0 36964.5 36995.0
[100] 37025.5 37056.0 37086.5 37117.5 37148.0 37178.5 37209.0 37239.5 37270.5
[109] 37300.0 37329.5 37360.0 37390.5 37421.0 37451.5 37482.5 37513.0 37543.5
[118] 37574.0 37604.5 37635.5 37665.0 37694.5 37725.0 37755.5 37786.0 37816.5
[127] 37847.5 37878.0 37908.5 37939.0 37969.5 38000.5 38030.5 38060.5 38091.0
[136] 38121.5 38152.0 38182.5 38213.5 38244.0 38274.5 38305.0 38335.5 38366.5
[145] 38396.0 38425.5 38456.0 38486.5 38517.0 38547.5 38578.5 38609.0 38639.5
[154] 38670.0 38700.5 38731.5 38761.0 38790.5 38821.0 38851.5 38882.0 38912.5
[163] 38943.5 38974.0 39004.5 39035.0 39065.5 39096.5 39126.0 39155.5 39186.0
[172] 39216.5 39247.0 39277.5 39308.5 39339.0 39369.5 39400.0 39430.5 39461.5
[181] 39491.5 39521.5 39552.0 39582.5 39613.0 39643.5 39674.5 39705.0 39735.5
[190] 39766.0 39796.5 39827.5 39857.0 39886.5 39917.0 39947.5 39978.0 40008.5
[199] 40039.5 40070.0 40100.5 40131.0 40161.5 40192.5 40222.0 40251.5 40282.0
[208] 40312.5 40343.0 40373.5 40404.5 40435.0 40465.5 40496.0 40526.5 40557.5
[217] 40587.0 40616.5 40647.0 40677.5 40708.0 40738.5 40769.5 40800.0 40830.5
[226] 40861.0 40891.5 40922.5 40952.5 40982.5 41013.0 41043.5 41074.0 41104.5
[235] 41135.5 41166.0 41196.5 41227.0 41257.5 41288.5 41318.0 41347.5 41378.0
[244] 41408.5 41439.0 41469.5 41500.5 41531.0 41561.5 41592.0 41622.5 41653.5
[253] 41683.0 41712.5 41743.0 41773.5 41804.0 41834.5 41865.5 41896.0 41926.5
[262] 41957.0 41987.5 42018.5 42048.0 42077.5 42108.0 42138.5 42169.0 42199.5
[271] 42230.5 42261.0 42291.5 42322.0 42352.5 42383.5 42413.5 42443.5 42474.0
[280] 42504.5 42535.0 42565.5 42596.5 42627.0 42657.5 42688.0 42718.5 42749.5
[289] 42779.0 42808.5 42839.0 42869.5 42900.0 42930.5 42961.5 42992.0 43022.5
[298] 43053.0 43083.5 43114.5 43144.0 43173.5 43204.0 43234.5 43265.0 43295.5
[307] 43326.5 43357.0 43387.5 43418.0 43448.5 43479.5 43509.0 43538.5 43569.0
[316] 43599.5 43630.0 43660.5 43691.5 43722.0 43752.5 43783.0 43813.5 43844.5
[325] 43874.5 43904.5 43935.0 43965.5 43996.0 44026.5 44057.5 44088.0 44118.5
[334] 44149.0 44179.5 44210.5 44240.0 44269.5 44300.0 44330.5 44361.0 44391.5
[343] 44422.5 44453.0 44483.5 44514.0 44544.5
tunits <- ncatt_get(ncin,"time","units")
nt <- dim(time)
nt
[1] 347
tunits
$hasatt
[1] TRUE

$value
[1] "days since 1900-01-01 00:00:00"
# Get temperature
dname <- "bottomT"

temp_array <- ncvar_get(ncin,dname)
dlname <- ncatt_get(ncin,dname,"long_name")
dunits <- ncatt_get(ncin,dname,"units")
fillvalue <- ncatt_get(ncin,dname,"_FillValue")
dim(temp_array)
[1] 519 498 347
# Get global attributes
title <- ncatt_get(ncin,0,"title")
institution <- ncatt_get(ncin,0,"institution")
datasource <- ncatt_get(ncin,0,"source")
references <- ncatt_get(ncin,0,"references")
history <- ncatt_get(ncin,0,"history")
Conventions <- ncatt_get(ncin,0,"Conventions")

# Convert time: split the time units string into fields
tustr <- strsplit(tunits$value, " ")
tdstr <- strsplit(unlist(tustr)[3], "-")
tmonth <- as.integer(unlist(tdstr)[2])
tday <- as.integer(unlist(tdstr)[3])
tyear <- as.integer(unlist(tdstr)[1])

# Here I deviate from the guide a little bit. Save this info:
dates <- chron(time, origin = c(tmonth, tday, tyear))

# Crop the date variable
months <- as.numeric(substr(dates, 2, 3))
years <- as.numeric(substr(dates, 8, 9))
years <- ifelse(years > 90, 1900 + years, 2000 + years)

# Replace netCDF fill values with NA's
temp_array[temp_array == fillvalue$value] <- NA

# Loop through all "dates", put into a list 
dlist <- list()

for(i in 1:length(months)) {
  
  temp_sub <- temp_array[, , i]
  
  dlist[[i]] <- temp_sub
  
}

# Name the list
names(dlist) <- paste(months, years, sep = "_")
str(dlist)
List of 347
 $ 2_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_1993: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_1993: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_1993: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_1994: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_1994: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_1994: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_1995: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_1995: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_1995: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_1996: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_1996: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_1996: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_1997: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_1997: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_1997: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_1998: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_1998: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_1998: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_1999: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_1999: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_1999: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_2000: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_2000: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_2000: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
  [list output truncated]
# Create data holding object
temp_data_list <- list()

# Loop through each month_year and extract raster values for the cpue data points
for(i in unique(dat$month_year)) { # We can use q1 as looping index, doesn't matter!
  
  # Set plot limits
  ymin = 54; ymax = 58; xmin = 12; xmax = 22
  
  # Subset a month-year combination
  temp_slice <- dlist[[i]]
  
  # Create raster for that year (i)
  r <- raster(t(temp_slice), xmn = min(lon), xmx = max(lon), ymn = min(lat), ymx = max(lat),
              crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
  
  # Flip...
  r <- flip(r, direction = 'y')
  
  plot(r, main = paste(i), ylim = c(ymin, ymax), xlim = c(xmin, xmax))
  
  # Filter the same year (i) in the cpue data and select only coordinates
  d_slice <- dat %>% filter(month_year == i) %>% dplyr::select(lon, lat)
  
  # Make into a SpatialPoints object
  data_sp <- SpatialPoints(d_slice)
  
  # Extract raster value (oxygen)
  rasValue <- raster::extract(r, data_sp)
  
  # Now we want to plot the results of the raster extractions by plotting the cpue data points over a raster and saving it for each year.
  # Make the SpatialPoints object into a raster again (for plot)
  df <- as.data.frame(data_sp)
  
  # Add in the raster value in the df holding the coordinates for the cpue data
  d_slice$temp <- rasValue
  
  # Add in which year
  d_slice$month_year <- i
  
  # Add each years' data in the list
  temp_data_list[[i]] <- d_slice
  
}
filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

# Now create a data frame from the list of all annual values
big_dat_temp <- dplyr::bind_rows(temp_data_list)

Salinity

# Open the netCDF file
ncin <- nc_open(paste0(home, "/data/NEMO/cmems_mod_bal_phy_my_P1M-m_1691065974147.nc"))

print(ncin)
File /Users/maxlindmark/Dropbox/Max work/R/pred-prey-overlap/data/NEMO/cmems_mod_bal_phy_my_P1M-m_1691065974147.nc (NC_FORMAT_CLASSIC):

     1 variables (excluding dimension variables):
        float sob[lon,lat,time]   
            standard_name: sea_water_salinity_at_sea_floor
            long_name: Sea water salinity at sea floor
            units: 0.001
            _FillValue: -999
            missing_value: -999
            cell_methods: time: mean
            unit_long: 0.001
            interval_write: 1 d
            interval_operation: 90 s
            online_operation: average
            _ChunkSizes: 1
             _ChunkSizes: 774
             _ChunkSizes: 763

     3 dimensions:
        time  Size:347 
            standard_name: time
            long_name: time
            bounds: time_bnds
            units: days since 1900-01-01 00:00:00
            calendar: standard
            axis: T
            _ChunkSizes: 512
            _CoordinateAxisType: Time
            valid_min: 34013
            valid_max: 44544.5
        lat  Size:498 
            standard_name: latitude
            long_name: Latitude of each location
            units: degrees_north
            axis: Y
            _CoordinateAxisType: Lat
            valid_min: 53.0082969665527
            valid_max: 61.2914962768555
        lon  Size:519 
            standard_name: longitude
            long_name: Longitude of each location
            units: degrees_east
            axis: X
            _CoordinateAxisType: Lon
            valid_min: 10.5414848327637
            valid_max: 24.9298248291016

    22 global attributes:
        CDI: Climate Data Interface version 1.9.9rc1 (https://mpimet.mpg.de/cdi)
        Conventions: CF-1.0
        source: CMEMS BAL MFC NEMO model output converted to NetCDF
        institution: Baltic MFC, PU Danish Meteorological Institute
        comment: Data on cropped native product grid. Horizontal velocities destagged
        grid_resolution: ~1 nautical mile (1min latitude; 1min40sec longitude)
        contact: servicedesk.cmems@mercator-ocean.eu
        references: https://marine.copernicus.eu/
        file_quality_index: 1
        compression: yes
        stop_date: 2021-12-01 12:00:00
        creation_date: 2023-03-14 13:30:28
        FROM_ORIGINAL_FILE__netcdf_version_id: 4.6.1 of Oct 19 2018 20:03:53 $
        FROM_ORIGINAL_FILE__easternmost_longitude: 30.2074012756348
        FROM_ORIGINAL_FILE__northernmost_latitude: 65.8909912109375
        FROM_ORIGINAL_FILE__westernmost_longitude: 9.04154205322266
        FROM_ORIGINAL_FILE__southernmost_latitude: 53.0082969665527
        start_date: 2021-12-01 12:00:00
        CDO: Climate Data Operators version 1.9.9rc1 (https://mpimet.mpg.de/cdo)
        title: CMEMS NEMO monthly integrated model fields
        _CoordSysBuilder: ucar.nc2.dataset.conv.CF1Convention
        history: Data extracted from dataset http://localhost:8080/thredds/dodsC/cmems_mod_bal_phy_my_P1M-m
# Get longitude and latitude
lon <- ncvar_get(ncin,"lon")
nlon <- dim(lon)
head(lon)
[1] 10.54148 10.56926 10.59704 10.62481 10.65259 10.68037
lat <- ncvar_get(ncin,"lat")
nlat <- dim(lat)
head(lat)
[1] 53.00830 53.02496 53.04163 53.05830 53.07496 53.09163
# Get time
time <- ncvar_get(ncin,"time")
time
  [1] 34013.0 34042.5 34073.0 34103.5 34134.0 34164.5 34195.5 34226.0 34256.5
 [10] 34287.0 34317.5 34348.5 34378.0 34407.5 34438.0 34468.5 34499.0 34529.5
 [19] 34560.5 34591.0 34621.5 34652.0 34682.5 34713.5 34743.0 34772.5 34803.0
 [28] 34833.5 34864.0 34894.5 34925.5 34956.0 34986.5 35017.0 35047.5 35078.5
 [37] 35108.5 35138.5 35169.0 35199.5 35230.0 35260.5 35291.5 35322.0 35352.5
 [46] 35383.0 35413.5 35444.5 35474.0 35503.5 35534.0 35564.5 35595.0 35625.5
 [55] 35656.5 35687.0 35717.5 35748.0 35778.5 35809.5 35839.0 35868.5 35899.0
 [64] 35929.5 35960.0 35990.5 36021.5 36052.0 36082.5 36113.0 36143.5 36174.5
 [73] 36204.0 36233.5 36264.0 36294.5 36325.0 36355.5 36386.5 36417.0 36447.5
 [82] 36478.0 36508.5 36539.5 36569.5 36599.5 36630.0 36660.5 36691.0 36721.5
 [91] 36752.5 36783.0 36813.5 36844.0 36874.5 36905.5 36935.0 36964.5 36995.0
[100] 37025.5 37056.0 37086.5 37117.5 37148.0 37178.5 37209.0 37239.5 37270.5
[109] 37300.0 37329.5 37360.0 37390.5 37421.0 37451.5 37482.5 37513.0 37543.5
[118] 37574.0 37604.5 37635.5 37665.0 37694.5 37725.0 37755.5 37786.0 37816.5
[127] 37847.5 37878.0 37908.5 37939.0 37969.5 38000.5 38030.5 38060.5 38091.0
[136] 38121.5 38152.0 38182.5 38213.5 38244.0 38274.5 38305.0 38335.5 38366.5
[145] 38396.0 38425.5 38456.0 38486.5 38517.0 38547.5 38578.5 38609.0 38639.5
[154] 38670.0 38700.5 38731.5 38761.0 38790.5 38821.0 38851.5 38882.0 38912.5
[163] 38943.5 38974.0 39004.5 39035.0 39065.5 39096.5 39126.0 39155.5 39186.0
[172] 39216.5 39247.0 39277.5 39308.5 39339.0 39369.5 39400.0 39430.5 39461.5
[181] 39491.5 39521.5 39552.0 39582.5 39613.0 39643.5 39674.5 39705.0 39735.5
[190] 39766.0 39796.5 39827.5 39857.0 39886.5 39917.0 39947.5 39978.0 40008.5
[199] 40039.5 40070.0 40100.5 40131.0 40161.5 40192.5 40222.0 40251.5 40282.0
[208] 40312.5 40343.0 40373.5 40404.5 40435.0 40465.5 40496.0 40526.5 40557.5
[217] 40587.0 40616.5 40647.0 40677.5 40708.0 40738.5 40769.5 40800.0 40830.5
[226] 40861.0 40891.5 40922.5 40952.5 40982.5 41013.0 41043.5 41074.0 41104.5
[235] 41135.5 41166.0 41196.5 41227.0 41257.5 41288.5 41318.0 41347.5 41378.0
[244] 41408.5 41439.0 41469.5 41500.5 41531.0 41561.5 41592.0 41622.5 41653.5
[253] 41683.0 41712.5 41743.0 41773.5 41804.0 41834.5 41865.5 41896.0 41926.5
[262] 41957.0 41987.5 42018.5 42048.0 42077.5 42108.0 42138.5 42169.0 42199.5
[271] 42230.5 42261.0 42291.5 42322.0 42352.5 42383.5 42413.5 42443.5 42474.0
[280] 42504.5 42535.0 42565.5 42596.5 42627.0 42657.5 42688.0 42718.5 42749.5
[289] 42779.0 42808.5 42839.0 42869.5 42900.0 42930.5 42961.5 42992.0 43022.5
[298] 43053.0 43083.5 43114.5 43144.0 43173.5 43204.0 43234.5 43265.0 43295.5
[307] 43326.5 43357.0 43387.5 43418.0 43448.5 43479.5 43509.0 43538.5 43569.0
[316] 43599.5 43630.0 43660.5 43691.5 43722.0 43752.5 43783.0 43813.5 43844.5
[325] 43874.5 43904.5 43935.0 43965.5 43996.0 44026.5 44057.5 44088.0 44118.5
[334] 44149.0 44179.5 44210.5 44240.0 44269.5 44300.0 44330.5 44361.0 44391.5
[343] 44422.5 44453.0 44483.5 44514.0 44544.5
tunits <- ncatt_get(ncin,"time","units")
nt <- dim(time)
nt
[1] 347
tunits
$hasatt
[1] TRUE

$value
[1] "days since 1900-01-01 00:00:00"
# Get Salinity
dname <- "sob"

sal_array <- ncvar_get(ncin,dname)
dlname <- ncatt_get(ncin,dname,"long_name")
dunits <- ncatt_get(ncin,dname,"units")
fillvalue <- ncatt_get(ncin,dname,"_FillValue")
dim(sal_array)
[1] 519 498 347
# Get global attributes
title <- ncatt_get(ncin,0,"title")
institution <- ncatt_get(ncin,0,"institution")
datasource <- ncatt_get(ncin,0,"source")
references <- ncatt_get(ncin,0,"references")
history <- ncatt_get(ncin,0,"history")
Conventions <- ncatt_get(ncin,0,"Conventions")

# Convert time: split the time units string into fields
tustr <- strsplit(tunits$value, " ")
tdstr <- strsplit(unlist(tustr)[3], "-")
tmonth <- as.integer(unlist(tdstr)[2])
tday <- as.integer(unlist(tdstr)[3])
tyear <- as.integer(unlist(tdstr)[1])

# Here I deviate from the guide a little bit. Save this info:
dates <- chron(time, origin = c(tmonth, tday, tyear))

# Crop the date variable
months <- as.numeric(substr(dates, 2, 3))
years <- as.numeric(substr(dates, 8, 9))
years <- ifelse(years > 90, 1900 + years, 2000 + years)

# Replace netCDF fill values with NA's
sal_array[sal_array == fillvalue$value] <- NA

# Loop through all "dates", put into a list 
dlist <- list()

for(i in 1:length(months)) {
  
  sal_sub <- sal_array[, , i]
  
  dlist[[i]] <- sal_sub
  
}

# Name the list
names(dlist) <- paste(months, years, sep = "_")
str(dlist)
List of 347
 $ 2_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_1993: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_1993: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_1993: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_1994: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_1994: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_1994: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_1995: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_1995: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_1995: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_1996: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_1996: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_1996: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_1997: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_1997: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_1997: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_1998: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_1998: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_1998: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_1999: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_1999: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_1999: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 5_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 6_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 7_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 8_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 9_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 10_2000: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 11_2000: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 12_2000: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
  [list output truncated]
# Create data holding object
sal_data_list <- list()

# Loop through each month_year and extract raster values for the cpue data points
for(i in unique(dat$month_year)) { # We can use q1 as looping index, doesn't matter!
  
  # Set plot limits
  ymin = 54; ymax = 58; xmin = 12; xmax = 22
  
  # Subset a month-year combination
  sal_slice <- dlist[[i]]
  
  # Create raster for that year (i)
  r <- raster(t(sal_slice), xmn = min(lon), xmx = max(lon), ymn = min(lat), ymx = max(lat),
              crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
  
  # Flip...
  r <- flip(r, direction = 'y')
  
  plot(r, main = paste(i), ylim = c(ymin, ymax), xlim = c(xmin, xmax))
  
  # Filter the same year (i) in the cpue data and select only coordinates
  d_slice <- dat %>% filter(month_year == i) %>% dplyr::select(lon, lat)
  
  # Make into a SpatialPoints object
  data_sp <- SpatialPoints(d_slice)
  
  # Extract raster value (oxygen)
  rasValue <- raster::extract(r, data_sp)
  
  # Now we want to plot the results of the raster extractions by plotting the cpue data points over a raster and saving it for each year.
  # Make the SpatialPoints object into a raster again (for plot)
  df <- as.data.frame(data_sp)
  
  # Add in the raster value in the df holding the coordinates for the cpue data
  d_slice$sal <- rasValue
  
  # Add in which year
  d_slice$month_year <- i
  
  # Add each years' data in the list
  sal_data_list[[i]] <- d_slice
  
}
filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

filter: removed 931,209 rows (98%), 16,337 rows remaining

# Now create a data frame from the list of all annual values
big_dat_sal <- dplyr::bind_rows(sal_data_list)
env_dat <- left_join(big_dat_oxy, big_dat_temp, by = c("month_year", "lon", "lat")) |> 
  left_join(big_dat_sal, by = c("month_year", "lon", "lat"))
left_join: added one column (temp)
           > rows only in x         0
           > rows only in y  (      0)
           > matched rows     947,546
           >                 =========
           > rows total       947,546
left_join: added one column (sal)
           > rows only in x         0
           > rows only in y  (      0)
           > matched rows     947,546
           >                 =========
           > rows total       947,546
# Now join these data with the full_dat
dat_full <- left_join(dat, env_dat, by = c("month_year", "lon", "lat"))
left_join: added 3 columns (oxy, temp, sal)
           > rows only in x         0
           > rows only in y  (      0)
           > matched rows     947,546
           >                 =========
           > rows total       947,546

Add ICES areas

# https://stackoverflow.com/questions/34272309/extract-shapefile-value-to-point-with-r
# https://gis.ices.dk/sf/
shape <- shapefile(paste0(home, "/data/ICES-StatRec-mapto-ICES-Areas/StatRec_map_Areas_Full_20170124.shp"))
head(shape)
  OBJECTID ID ICESNAME SOUTH WEST NORTH EAST AREA_KM2 stat_x stat_y Area_27
1        1 47     47A0  59.0  -44  59.5  -43     3178  -43.5  59.25  14.b.2
2        2 48     48A0  59.5  -44  60.0  -43     3132  -43.5  59.75  14.b.2
3        3 49     49A0  60.0  -44  60.5  -43     3085  -43.5  60.25  14.b.2
4        4 50     50A0  60.5  -44  61.0  -43     3038  -43.5  60.75  14.b.2
5        5 51     51A0  61.0  -44  61.5  -43     2991  -43.5  61.25  14.b.2
6        6 52     52A0  61.5  -44  62.0  -43     2944  -43.5  61.75  14.b.2
          Perc      MaxPer RNDMaxPer AreasList Shape_Leng Shape_Area
1 100.00000000 100.0000000       100    14.b.2          3        0.5
2  84.12674145  84.1267414        84    14.b.2          3        0.5
3  24.99803694  24.9980369        25    14.b.2          3        0.5
4  11.97744244  11.9774424        12    14.b.2          3        0.5
5   3.89717625   3.8971762         4    14.b.2          3        0.5
6   0.28578428   0.2857843         0    14.b.2          3        0.5
pts <- SpatialPoints(cbind(dat_full$lon, dat_full$lat), 
                     proj4string = CRS(proj4string(shape)))

dat_full$subdiv <- over(pts, shape)$Area_27

# Rename subdivisions to the more common names and do some more filtering (by sub div and area)
sort(unique(dat_full$subdiv))
[1] "3.d.24"   "3.d.25"   "3.d.26"   "3.d.27"   "3.d.28.2"
dat_full <- dat_full |> 
  mutate(sub_div = factor(subdiv),
         sub_div = fct_recode(subdiv,
                              "24" = "3.d.24",
                              "25" = "3.d.25",
                              "26" = "3.d.26",
                              "27" = "3.d.27",
                              "28" = "3.d.28.1",
                              "28" = "3.d.28.2",
                              "29" = "3.d.29"),
         sub_div = as.character(sub_div)) |> 
  filter(sub_div %in% c("24", "25", "26", "27", "28", 2)) |> 
  filter(lat > 54 & lat < 59 & lon < 22)
Warning: There was 1 warning in `.fun()`.
ℹ In argument: `sub_div = fct_recode(...)`.
Caused by warning:
! Unknown levels in `f`: 3.d.28.1, 3.d.29
mutate: new variable 'sub_div' (character) with 5 unique values and 0% NA
filter: no rows removed
filter: no rows removed
# Add ICES rectangles
dat_full$ices_rect <- mapplots::ices.rect2(lon = dat_full$lon, lat = dat_full$lat)

plot_map +
  geom_raster(data = filter(dat_full, year == 1999), aes(X*1000, Y*1000, fill = oxy)) +
  facet_wrap(~sub_div)
filter: removed 914,872 rows (97%), 32,674 rows remaining
Warning: Removed 1524 rows containing missing values (`geom_raster()`).

dat_full <- dat_full |> dplyr::select(-subdiv)

Save

# Remove variables and save
pred_grid_93_06 <- dat_full |> filter(year < 2007)
filter: removed 490,110 rows (52%), 457,436 rows remaining
pred_grid_07_19 <- dat_full |> filter(year > 2006)
filter: removed 457,436 rows (48%), 490,110 rows remaining
write_csv(pred_grid_93_06, file = paste0(home, "/data/clean/pred_grid_(1_2).csv"))
write_csv(pred_grid_07_19, file = paste0(home, "/data/clean/pred_grid_(2_2).csv"))